{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Content under Creative Commons Attribution license CC-BY 4.0, code under BSD 3-Clause License © 2018 by D. Koehn, notebook style sheet by L.A. Barba, N.C. Clementi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<link href=\"https://fonts.googleapis.com/css?family=Merriweather:300,300i,400,400i,700,700i,900,900i\" rel='stylesheet' >\n",
       "<link href=\"https://fonts.googleapis.com/css?family=Source+Sans+Pro:300,300i,400,400i,700,700i\" rel='stylesheet' >\n",
       "<link href='http://fonts.googleapis.com/css?family=Source+Code+Pro:300,400' rel='stylesheet' >\n",
       "<style>\n",
       "\n",
       "@font-face {\n",
       "    font-family: \"Computer Modern\";\n",
       "    src: url('http://mirrors.ctan.org/fonts/cm-unicode/fonts/otf/cmunss.otf');\n",
       "}\n",
       "\n",
       "\n",
       "#notebook_panel { /* main background */\n",
       "    background: rgb(245,245,245);\n",
       "}\n",
       "\n",
       "div.cell { /* set cell width */\n",
       "    width: 800px;\n",
       "}\n",
       "\n",
       "div #notebook { /* centre the content */\n",
       "    background: #fff; /* white background for content */\n",
       "    width: 1000px;\n",
       "    margin: auto;\n",
       "    padding-left: 0em;\n",
       "}\n",
       "\n",
       "#notebook li { /* More space between bullet points */\n",
       "margin-top:0.5em;\n",
       "}\n",
       "\n",
       "/* draw border around running cells */\n",
       "div.cell.border-box-sizing.code_cell.running { \n",
       "    border: 1px solid #111;\n",
       "}\n",
       "\n",
       "/* Put a solid color box around each cell and its output, visually linking them*/\n",
       "div.cell.code_cell {\n",
       "    background-color: rgb(256,256,256); \n",
       "    border-radius: 0px; \n",
       "    padding: 0.5em;\n",
       "    margin-left:1em;\n",
       "    margin-top: 1em;\n",
       "}\n",
       "\n",
       "\n",
       "div.text_cell_render{\n",
       "    font-family: 'Source Sans Pro', sans-serif;\n",
       "    line-height: 140%;\n",
       "    font-size: 110%;\n",
       "    width:680px;\n",
       "    margin-left:auto;\n",
       "    margin-right:auto;\n",
       "}\n",
       "\n",
       "/* Formatting for header cells */\n",
       ".text_cell_render h1 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-style:regular;\n",
       "    font-weight: bold;    \n",
       "    font-size: 250%;\n",
       "    line-height: 100%;\n",
       "    color: #004065;\n",
       "    margin-bottom: 1em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\t\n",
       ".text_cell_render h2 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-weight: bold; \n",
       "    font-size: 180%;\n",
       "    line-height: 100%;\n",
       "    color: #0096d6;\n",
       "    margin-bottom: 0.5em;\n",
       "    margin-top: 0.5em;\n",
       "    display: block;\n",
       "}\t\n",
       "\n",
       ".text_cell_render h3 {\n",
       "    font-family: 'Merriweather', serif;\n",
       "\tfont-size: 150%;\n",
       "    margin-top:12px;\n",
       "    margin-bottom: 3px;\n",
       "    font-style: regular;\n",
       "    color: #008367;\n",
       "}\n",
       "\n",
       ".text_cell_render h4 {    /*Use this for captions*/\n",
       "    font-family: 'Merriweather', serif;\n",
       "    font-weight: 300; \n",
       "    font-size: 100%;\n",
       "    line-height: 120%;\n",
       "    text-align: left;\n",
       "    width:500px;\n",
       "    margin-top: 1em;\n",
       "    margin-bottom: 2em;\n",
       "    margin-left: 80pt;\n",
       "    font-style: regular;\n",
       "}\n",
       "\n",
       ".text_cell_render h5 {  /*Use this for small titles*/\n",
       "    font-family: 'Source Sans Pro', sans-serif;\n",
       "    font-weight: regular;\n",
       "    font-size: 130%;\n",
       "    color: #e31937;\n",
       "    font-style: italic;\n",
       "    margin-bottom: .5em;\n",
       "    margin-top: 1em;\n",
       "    display: block;\n",
       "}\n",
       "\n",
       ".text_cell_render h6 { /*use this for copyright note*/\n",
       "    font-family: 'Source Code Pro', sans-serif;\n",
       "    font-weight: 300;\n",
       "    font-size: 9pt;\n",
       "    line-height: 100%;\n",
       "    color: grey;\n",
       "    margin-bottom: 1px;\n",
       "    margin-top: 1px;\n",
       "}\n",
       "\n",
       "    .CodeMirror{\n",
       "            font-family: \"Source Code Pro\";\n",
       "\t\t\tfont-size: 90%;\n",
       "    }\n",
       "/*    .prompt{\n",
       "        display: None;\n",
       "    }*/\n",
       "\t\n",
       "    \n",
       "    .warning{\n",
       "        color: rgb( 240, 20, 20 )\n",
       "        }  \n",
       "</style>\n",
       "<script>\n",
       "    MathJax.Hub.Config({\n",
       "                        TeX: {\n",
       "                           extensions: [\"AMSmath.js\"], \n",
       "                           equationNumbers: { autoNumber: \"AMS\", useLabelIds: true}\n",
       "                           },\n",
       "                tex2jax: {\n",
       "                    inlineMath: [ ['$','$'], [\"\\\\(\",\"\\\\)\"] ],\n",
       "                    displayMath: [ ['$$','$$'], [\"\\\\[\",\"\\\\]\"] ]\n",
       "                },\n",
       "                displayAlign: 'center', // Change this to 'center' to center equations.\n",
       "                \"HTML-CSS\": {\n",
       "                    styles: {'.MathJax_Display': {\"margin\": 4}}\n",
       "                }\n",
       "        });\n",
       "</script>\n"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Execute this cell to load the notebook's style sheet, then ignore it\n",
    "from IPython.core.display import HTML\n",
    "css_file = '../style/custom.css'\n",
    "HTML(open(css_file, \"r\").read())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# From 2D acoustic to 2D SH finite difference modelling\n",
    "\n",
    "In the last [lesson](http://nbviewer.jupyter.org/github/daniel-koehn/Theory-of-seismic-waves-II/blob/master/06_2D_SH_Love_wave_modelling/1_2D_SH_FD_staggered.ipynb), we derived a 2nd order space/time finite difference solution on a staggered Cartesian grid to solve the 2D isotropic elastic SH problem. In the next step we will implement this FD scheme starting from our 2D acoustic code ..."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## FD solution of 2D isotropic SH problem\n",
    "\n",
    "The 2nd order space/time FD scheme on a staggered grid \n",
    "\n",
    "<img src=\"images/SG_SH-Cart.png\" width=\"50%\">\n",
    "\n",
    "using explicit time-stepping with the Leapfrog method to solve the 2D SH problem can be written as \n",
    "\n",
    "\\begin{align}\n",
    "\\scriptsize \n",
    "v_y(i,j,n+1/2) &\\scriptsize= \\frac{dt}{\\rho}\\biggl\\{\\biggl(\\frac{\\partial \\sigma_{yx}}{\\partial x}\\biggr)^c(i,j,n) + \\biggl(\\frac{\\partial \\sigma_{yz}}{\\partial z}\\biggl)^c(i,j,n) + f_y(i,j,n)\\biggr\\} + v_y(i,j,n-1/2), \\notag\\\\\n",
    "\\scriptsize\\sigma_{yx}(i+1/2,j,n+1) &\\scriptsize= dt\\; \\mu_x \\biggl(\\frac{\\partial v_{y}}{\\partial x}\\biggr)^c(i+1/2,j,n+1/2) + \\sigma_{yx}(i+1/2,j,n),\\notag\\\\\n",
    "\\scriptsize\\sigma_{yz}(i,j+1/2,n+1) &\\scriptsize= dt\\; \\mu_z \\biggl(\\frac{\\partial v_{y}}{\\partial z}\\biggr)^c(i,j+1/2,n+1/2) + \\sigma_{yz}(i,j+1/2,n), \\notag\\\\\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "with the spatial derivatives on the RHS approximated by\n",
    "\n",
    "\\begin{equation}\n",
    "\\scriptsize \\biggl(\\frac{\\partial \\sigma_{yx}}{\\partial x}\\biggr)^c(i,j) \\approx \\frac{\\sigma_{yx}(i+1/2,j) - \\sigma_{yx}(i-1/2,j)}{dx},\\; \n",
    "\\biggl(\\frac{\\partial \\sigma_{yz}}{\\partial z}\\biggr)^c(i,j) \\approx \\frac{\\sigma_{yz}(i,j+1/2) - \\sigma_{yz}(i,j-1/2)}{dz},\\notag\n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\scriptsize \\biggl(\\frac{\\partial v_{y}}{\\partial x}\\biggr)^c(i+1/2,j) \\approx \\frac{v_y(i+1,j) - v_{y}(i,j)}{dx},\\notag \n",
    "\\end{equation}\n",
    "\n",
    "\\begin{equation}\n",
    "\\scriptsize \\biggl(\\frac{\\partial v_{y}}{\\partial z}\\biggr)^c(i,j+1/2) \\approx \\scriptsize \\frac{v_y(i,j+1) - v_{y}(i,j)}{dz}.\\notag \n",
    "\\end{equation}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "and the harmonically averaged shear moduli:\n",
    "\n",
    "\\begin{align}\n",
    "\\mu_{x} &= \\mu(i+1/2,j) = 2 \\biggl[\\biggl(\\mu^{-1}_{i,j}+\\mu^{-1}_{i+1,j}\\biggr)\\biggr]^{-1} \\notag\\\\\n",
    "\\mu_{z} &= \\mu(i,j+1/2) = 2 \\biggl[\\biggl(\\mu^{-1}_{i,j}+\\mu^{-1}_{i,j+1}\\biggr)\\biggr]^{-1} \\notag\n",
    "\\end{align}"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Initial and boundary conditions\n",
    "\n",
    "Because we have analytical solutions for wave propagation in homogeneous media, we should test our first code implementation for a similar medium, by setting density $\\rho$ and shear modulus $\\mu$ to constant values $\\rho_0,\\; \\mu_0$\n",
    "\n",
    "\\begin{align}\n",
    "\\rho(i,j) &= \\rho_0 \\notag \\\\\n",
    "\\mu(i,j) &= \\mu_0 = \\rho_0 V_{s0}^2\\notag\n",
    "\\end{align}\n",
    "\n",
    "at each spatial grid point $i = 0, 1, 2, ..., nx$; $j = 0, 1, 2, ..., nz$, in order to compare the numerical with the analytical solution. For a complete description of the problem we also have to define initial and boundary conditions. The **initial condition** is \n",
    "\n",
    "\\begin{equation}\n",
    "v_y(i,j,0) = \\sigma_{yx}(i+1/2,j,0) = \\sigma_{yz}(i,j+1/2,0) = 0, \\nonumber\n",
    "\\end{equation}\n",
    "\n",
    "so the modelling starts with zero particel velocity and shear stress amplitudes at each spatial grid point. As **boundary conditions**, we assume \n",
    "\n",
    "\\begin{align}\n",
    "v_y(0,j,n) &= \\sigma_{yx}(1/2,j,n) = \\sigma_{yz}(0,j+1/2,n) = 0, \\nonumber\\\\\n",
    "v_y(nx,j,n) &= \\sigma_{yx}(nx+1/2,j,n) = \\sigma_{yz}(nx,j+1/2,n) = 0, \\nonumber\\\\\n",
    "v_y(i,0,n) &= \\sigma_{yx}(i+1/2,0,n) = \\sigma_{yz}(i,1/2,n) = 0, \\nonumber\\\\\n",
    "v_y(i,nz,n) &= \\sigma_{yx}(i+1/2,nz,n) = \\sigma_{yz}(i,nz+1/2,n) = 0, \\nonumber\\\\\n",
    "\\end{align}\n",
    "\n",
    "for all time steps n. This **Dirichlet boundary condition**, leads to artifical boundary reflections which would obviously not describe a homogeneous medium. For now, we simply extend the model, so that boundary reflections are not recorded at the receiver positions.\n",
    "\n",
    "Let's implement it ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "code_folding": [
     0
    ]
   },
   "outputs": [],
   "source": [
    "# Import Libraries \n",
    "# ----------------\n",
    "import numpy as np\n",
    "from numba import jit\n",
    "import matplotlib\n",
    "import matplotlib.pyplot as plt\n",
    "from pylab import rcParams\n",
    "\n",
    "# Ignore Warning Messages\n",
    "# -----------------------\n",
    "import warnings\n",
    "warnings.filterwarnings(\"ignore\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "code_folding": [],
    "scrolled": false
   },
   "outputs": [],
   "source": [
    "# Definition of modelling parameters\n",
    "# ----------------------------------\n",
    "xmax = 500.0 # maximum spatial extension of the 1D model in x-direction (m)\n",
    "zmax = xmax  # maximum spatial extension of the 1D model in z-direction(m)\n",
    "#dx   = 1.0   # grid point distance in x-direction\n",
    "#dz   = dx    # grid point distance in z-direction\n",
    "\n",
    "tmax = 0.502   # maximum recording time of the seismogram (s)\n",
    "#dt   = 0.0010  # time step\n",
    "\n",
    "vs0  = 580.   # S-wave speed in medium (m/s)\n",
    "rho0 = 1000.  # Density in medium (kg/m^3)\n",
    "\n",
    "# acquisition geometry\n",
    "xr = 330.0 # x-receiver position (m)\n",
    "zr = xr    # z-receiver position (m)\n",
    "\n",
    "xsrc = 250.0 # x-source position (m)\n",
    "zsrc = 250.0 # z-source position (m)\n",
    "\n",
    "f0   = 40. # dominant frequency of the source (Hz)\n",
    "t0   = 4. / f0 # source time shift (s)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Comparison of 2D finite difference with analytical solution for homogeneous Vs model\n",
    "\n",
    "In a previous [exercise](https://danielkoehnsite.files.wordpress.com/2018/04/ex_02_tew21.pdf) you proved that the analytical solutions for the homogeneous 2D acoustic and 2D SH problem, beside a density factor $1/\\rho_0$, are actual identical . In the function below we solve the homogeneous 2D SH problem by centered 2nd order spatial/temporal difference operators and compare the numerical results with the analytical solution: \n",
    "\n",
    "\\begin{equation}\n",
    "u_{y,analy}(x,z,t) = G_{2D} * S \\nonumber \n",
    "\\end{equation}\n",
    "\n",
    "with the 2D Green's function:\n",
    "\n",
    "\\begin{equation}\n",
    "G_{2D}(x,z,t) = \\dfrac{1}{2\\pi \\rho_0 V_{s0}^2}\\dfrac{H\\biggl((t-t_s)-\\dfrac{|r|}{V_{s0}}\\biggr)}{\\sqrt{(t-t_s)^2-\\dfrac{r^2}{V_{s0}^2}}}, \\nonumber \n",
    "\\end{equation}\n",
    "\n",
    "where $H$ denotes the Heaviside function, $r = \\sqrt{(x-x_s)^2+(z-z_s)^2}$ the source-receiver distance (offset) and $S$ the source wavelet. Keep in mind that the stress-velocity code computes the **particle velocities** $\\mathbf{v_{y,analy}}$, while the analytical solution is expressed in terms of the **displacement** $\\mathbf{u_{y,analy}}$. Therefore, we have to take the first derivative of the analytical solution, before comparing the numerical with the analytical solution:\n",
    "\n",
    "\\begin{equation}\n",
    "v_{y,analy}(x,z,t) = \\frac{\\partial u_{y,analy}}{\\partial t} \\nonumber \n",
    "\\end{equation}\n",
    "\n",
    "To implement the 2D SH code, we first introduce functions to update the particle velocity $v_y$ ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Particle velocity vy update\n",
    "# ---------------------------\n",
    "@jit(nopython=True) # use JIT for C-performance\n",
    "def update_vel(vy, syx, syz, dx, dz, dt, nx, nz, rho):\n",
    "    \n",
    "    for i in range(1, nx - 1):\n",
    "        for j in range(1, nz - 1):\n",
    "            \n",
    "            # Calculate spatial derivatives            \n",
    "            syx_x = (syx[i,j] - syx[i - 1,j]) / dx\n",
    "            syz_z = (syz[i,j] - syz[i,j - 1]) / dz\n",
    "            \n",
    "            # Update particle velocities\n",
    "            vy[i,j] = vy[i,j] + (dt/rho[i,j]) * (syx_x + syz_z)                \n",
    "                    \n",
    "    return vy"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... update the shear stress components $\\sigma_{yx}$ and $\\sigma_{yz}$ ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Shear stress syx, syz updates\n",
    "# -----------------------------\n",
    "@jit(nopython=True) # use JIT for C-performance\n",
    "def update_stress(vy, syx, syz, dx, dz, dt, nx, nz, mux, muz):\n",
    "    \n",
    "    for i in range(1, nx - 1):\n",
    "        for j in range(1, nz - 1):\n",
    "            \n",
    "            # Calculate spatial derivatives            \n",
    "            vy_x = (vy[i + 1,j] - vy[i,j]) / dx\n",
    "            vy_z = (vy[i,j + 1] - vy[i,j]) / dz\n",
    "            \n",
    "            # Update shear stresses\n",
    "            syx[i,j] = syx[i,j] + dt * mux[i,j] * vy_x\n",
    "            syz[i,j] = syz[i,j] + dt * muz[i,j] * vy_z\n",
    "                    \n",
    "    return syx, syz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "... and harmonically averaging the shear modulus ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Harmonic averages of shear modulus\n",
    "# ----------------------------------\n",
    "@jit(nopython=True) # use JIT for C-performance\n",
    "def shear_avg(mu, nx, nz, mux, muz):\n",
    "    \n",
    "    for i in range(1, nx - 1):\n",
    "        for j in range(1, nz - 1):\n",
    "            \n",
    "            # Calculate harmonic averages of shear moduli        \n",
    "            mux[i,j] = 2 / (1 / mu[i + 1,j] + 1 / mu[i,j])\n",
    "            muz[i,j] = 2 / (1 / mu[i,j + 1] + 1 / mu[i,j])\n",
    "            \n",
    "    return mux, muz"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Finally, we can modify the main FD code ..."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "code_folding": [
     42
    ]
   },
   "outputs": [],
   "source": [
    "# 2D SH Wave Propagation (Finite Difference Solution) \n",
    "# ---------------------------------------------------\n",
    "def FD_2D_SH_JIT(dt,dx,dz,f0,xsrc,zsrc):\n",
    "        \n",
    "    nx = (int)(xmax/dx) # number of grid points in x-direction\n",
    "    print('nx = ',nx)\n",
    "    \n",
    "    nz = (int)(zmax/dz) # number of grid points in x-direction\n",
    "    print('nz = ',nz)\n",
    "            \n",
    "    nt = (int)(tmax/dt) # maximum number of time steps            \n",
    "    print('nt = ',nt)\n",
    "    \n",
    "    ir = (int)(xr/dx)      # receiver location in grid in x-direction    \n",
    "    jr = (int)(zr/dz)      # receiver location in grid in z-direction\n",
    "    \n",
    "    isrc = (int)(xsrc/dx)  # source location in grid in x-direction\n",
    "    jsrc = (int)(zsrc/dz)  # source location in grid in x-direction\n",
    "\n",
    "    # Source time function (Gaussian)\n",
    "    # -------------------------------\n",
    "    src  = np.zeros(nt + 1)\n",
    "    time = np.linspace(0 * dt, nt * dt, nt)\n",
    "\n",
    "    # 1st derivative of a Gaussian\n",
    "    src  = -2. * (time - t0) * (f0 ** 2) * (np.exp(- (f0 ** 2) * (time - t0) ** 2))\n",
    "\n",
    "    # Analytical solution\n",
    "    # -------------------\n",
    "    G        = time * 0.\n",
    "    vy_analy = time * 0.\n",
    "\n",
    "    # Initialize coordinates\n",
    "    # ----------------------\n",
    "    x    = np.arange(nx)\n",
    "    x    = x * dx       # coordinates in x-direction (m)\n",
    "\n",
    "    z    = np.arange(nz)\n",
    "    z    = z * dz       # coordinates in z-direction (m)\n",
    "    \n",
    "    # calculate source-receiver distance\n",
    "    r = np.sqrt((x[ir] - x[isrc])**2 + (z[jr] - z[jsrc])**2)\n",
    "    \n",
    "    for it in range(nt): # Calculate Green's function (Heaviside function)\n",
    "        if (time[it] - r / vs0) >= 0:\n",
    "            G[it] = 1. / (2 * np.pi * rho0 * vs0**2) * (1. / np.sqrt(time[it]**2 - (r/vs0)**2))\n",
    "    Gc   = np.convolve(G, src * dt)\n",
    "    Gc   = Gc[0:nt]\n",
    "    \n",
    "    # compute vy_analy from uy_analy\n",
    "    for i in range(1, nt - 1):\n",
    "        vy_analy[i] = (Gc[i+1] - Gc[i-1]) / (2.0 * dt)           \n",
    "    \n",
    "    # Initialize empty pressure arrays\n",
    "    # --------------------------------\n",
    "    vy    = np.zeros((nx,nz)) # particle velocity vy\n",
    "    syx   = np.zeros((nx,nz)) # shear stress syx\n",
    "    syz   = np.zeros((nx,nz)) # shear stress syz        \n",
    "\n",
    "    # Initialize model (assume homogeneous model)\n",
    "    # -------------------------------------------    \n",
    "    vs    = np.zeros((nx,nz))\n",
    "    vs    = vs + vs0      # initialize wave velocity in model\n",
    "    \n",
    "    rho   = np.zeros((nx,nz))\n",
    "    rho   = rho + rho0    # initialize wave velocity in model\n",
    "    \n",
    "    # calculate shear modulus\n",
    "    mu    = np.zeros((nx,nz))\n",
    "    mu    = rho * vs ** 2 \n",
    "    \n",
    "    # harmonic average of shear moduli\n",
    "    # --------------------------------\n",
    "    mux   = mu # initialize harmonic average mux \n",
    "    muz   = mu # initialize harmonic average muz\n",
    "\n",
    "    mux, muz = shear_avg(mu, nx, nz, mux, muz)\n",
    "    \n",
    "    # Initialize empty seismogram\n",
    "    # ---------------------------\n",
    "    seis = np.zeros(nt) \n",
    "    \n",
    "    # Time looping\n",
    "    # ------------\n",
    "    for it in range(nt):\n",
    "    \n",
    "        # Update particle velocity vy\n",
    "        # ---------------------------\n",
    "        vy = update_vel(vy, syx, syz, dx, dz, dt, nx, nz, rho)\n",
    "\n",
    "        # Add Source Term at (isrc,jsrc)\n",
    "        # ------------------------------\n",
    "        # Absolute particle velocity w.r.t analytical solution\n",
    "        vy[isrc,jsrc] = vy[isrc,jsrc] + (dt * src[it] / (rho[isrc,jsrc] * dx * dz))\n",
    "        \n",
    "        # Update shear stress syx, syz\n",
    "        # ----------------------------\n",
    "        syx, syz = update_stress(vy, syx, syz, dx, dz, dt, nx, nz, mux, muz)                \n",
    "                           \n",
    "        # Output of Seismogram\n",
    "        # -----------------\n",
    "        seis[it] = vy[ir,jr]   \n",
    "        \n",
    "    # Compare FD Seismogram with analytical solution\n",
    "    # ---------------------------------------------- \n",
    "    # Define figure size\n",
    "    rcParams['figure.figsize'] = 12, 5\n",
    "    plt.plot(time, seis, 'b-',lw=3,label=\"FD solution\") # plot FD seismogram\n",
    "    Analy_seis = plt.plot(time,vy_analy,'r--',lw=3,label=\"Analytical solution\") # plot analytical solution\n",
    "    plt.xlim(time[0], time[-1])\n",
    "    plt.title('Seismogram')\n",
    "    plt.xlabel('Time (s)')\n",
    "    plt.ylabel('Amplitude')\n",
    "    plt.legend()\n",
    "    plt.grid()\n",
    "    plt.show()     \n",
    "    \n",
    "    return mux, muz"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "nx =  500\n",
      "nz =  500\n",
      "nt =  502\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAtUAAAFNCAYAAADPbpPYAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDIuMi4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvhp/UCwAAIABJREFUeJzs3Xd81dX9x/HXudkLwgjzhg1hJBBEQAUREfeo4qgUF9piHR3aOn511NlaRx2l2rpKrYq4cC9QooKCCgRky07YK5C97vn9cS/3ZkEuZHyT3Pfz8bjme77fe873E3JMPjk533OMtRYRERERETl6LqcDEBERERFp7pRUi4iIiIjUkZJqEREREZE6UlItIiIiIlJHSqpFREREROpISbWIiIiISB0pqRYRaSKMMXnGmF5OxyEiIkdOSbWISD0zxow2xnxjjNlvjNlrjJlnjBleWz1rbby1dn1jxCgiIvUr3OkARERaEmNMK+AD4DrgdSASOBEodjKu+mCMCbPWljsdh4hIU6SRahGR+tUPwFo73Vpbbq0ttNZ+Zq1dCmCMudoYs9IYs88Y86kxpvvBisYYa4zp4zs+yxizwhiTa4zZYoz5o+/8WGNMtjHmVmPMTmPMNmPM+b73r/GNjP+pQptRxpgnjDFbfa8njDFRFa7f6mtjqzHml1VimGaMecYY85ExJh842RhztjFmsTHmgDEmyxhzT4W2evjqT/Zd22eM+bUxZrgxZqkxJscYM7Vh//lFRJyhpFpEpH6tAcqNMf81xpxpjGlz8IIx5nzgT8AEIAn4Gph+iHZeAK611iYAqcAXFa51AqKBrsDdwHPAZcAwvKPid1eYm30HcByQDgwBRgB3+uI5A7gZGA/0AU6qIY5fAA8CCcBcIB+4AkgEzgau831eFY0E+gI/B57wxTAeGARcYoyp6T4iIs1as0uqjTEv+kZnltVDWycbYzIrvIpq+OEgIhI0a+0BYDRg8Sa7u4wx7xljOgLXAn+11q601pYBfwHSK45WV1AKDDTGtLLW7rPWLqpy7UFrbSnwGtAeeNJam2utXQ4sBwb73jsJuM9au9Nauwu4F7jcd+0S4D/W2uXW2gLftaretdbOs9Z6rLVF1toMa+2PvvJSvL8UVE2S7/e99zO8Sfh03/234P1FYmhw/5oiIs1Hs0uqgWnAGfXRkLV2jrU23VqbDowDCoDP6qNtEQldvqT5KmutG+8ocxe8I7bdgSd90yBygL2AwTviXNWFwFnAJmPMl8aY4ytc21NhbnOh7+OOCtcLgXjfcRdgU4Vrm3znDl7LqnCt4nGN54wxI40xc4wxu4wx+4Ff403qK6oay6FiExFpMZpdUm2t/QrvDyI/Y0xvY8wnxpiFxpivjTH9j6Lpi4CPfaM1IiL1wlq7Cu9gQCreBPVaa21ihVeMtfabGup9b639GdABeAfvQ49HYyveZP6gbr5zANsAd4VryTV9ClXKrwLvAcnW2tbAv/D+YiAiEtKaXVJ9CM8Cv7HWDgP+CDx9FG1cyqHnNoqIBMUY098Y8wdjjNtXTgYmAvPxJqD/Z4wZ5LvW2hhzcQ1tRBpjJhljWvumeBwAjnbVjenAncaYJGNMe7xzsF/2XXsdmGyMGWCMifVdq00CsNdaW2SMGYF3zrWISMhr9kvqGWPigROAN4zxD5ZE+a5NAO6rodoWa+3pFdroDKQBnzZstCISAnLxPqh3szEmEcjBu8TeLdbaA77vWa/55lHvB2YBb9TQzuXAVGNMGLAa74OIR+MBoBWw1Fd+w3cOa+3HxpingDmAB7jfd9/DLf93PfCYbxWPL/Em5olHGZuISIthrK36l72mzxjTA/jAWpvqWxN2tbW2cx3a+x0wyFo7pZ5CFBFpdowxA4BlQJTvQUoREQlSs5/+4XvSfsPBP6EaryFH2MxENPVDREKQMeYC33STNsDfgPeVUIuIHLlml1QbY6YD3wIpvg0QrsG7ZNQ1xpgleJeS+tkRtNcD78M5X9Z/tCIiTd61wC5gHd5529c5G46ISPPULKd/iIiIiIg0Jc1upFpEREREpKlRUi0iIiIiUkfNakm9xMRE26dPH6fDkCYkPz+fuLg4p8OQJkb9QqpSn5Cq1CekJhX7xcKFC3dba5OCrduskuqOHTvyww8/OB2GNCEZGRmMHTvW6TCkiVG/kKrUJ6Qq9QmpScV+YYzZdCR1Nf1DRERERKSOlFSLiIiIiNSRkmoRERERkTpybE61MSYa+AqI8sXxprX2z07FIyIiIqGptLSU7OxsioqKnA5FHBAdHY3b7SYiIqJO7Tj5oGIxMM5am2eMiQDmGmM+ttbOdzAmERERCTHZ2dkkJCTQo0cPjDFOhyONyFrLnj17yM7OpmfPnnVqy7HpH9Yrz1eM8L20vaOIiIg0qqKiItq1a6eEOgQZY2jXrl29/JXC0TnVxpgwY0wmsBOYZa1d4GQ8IiIiEpqUUIeu+vraG2udHxw2xiQCM4HfWGuXVbk2BZgCkJSUNOz11193IEJpqvLy8oiPj3c6DGli1C+kKvUJqapin2jdujVOby6XmJjIoEGD/OVXX32VzZs3M3HiRHr06EFBQQEdOnTgd7/7HWeeeWad7rVp0yYuueQSFiw4/Fjmo48+yh//+Ed/efz48cyePbtO926q1q5dy/79+yv1i5NPPnmhtfbYYNtoEpu/WGtzjDEZwBnAsirXngWeBUhJSbFaqF0q0uL9UhP1i8ZlLezZAwkJEBXldDQ1U5+Qqir2iZUrV5KQkOBoPDExMSxdurTSud27d3PiiSfywQcfAJCZmcn5559Pu3btOOWUU476XvHx8bhcrlo/58cee4x7773XX64tCW/OoqOjGTp0aJ2+Vzg2/cMYk+QbocYYEwOMB1Y5FY+IiByZsjJ48EHo2hX6Ju2ja9tCrr4adu50OjKRlik9PZ27776bqVOnVrv25Zdfkp6eTnp6OkOHDiU3NxdrLbfccgupqamkpaUxY8aMavWmTZvGjTfe6C+fc845ZGRkcPvtt1NYWEh6ejqTJk0C8I/gHqrdgwnpRRddRP/+/Zk0aRJNYUZEY3FypLoz8F9jTBje5P51a+0HDsYjIiJBOnAAzjvX0vOracziUQaxguKCSD7/zylcMedxXpibQteuTkcp0nwcTGABevbsycyZM2t83zHHHMMjjzxS7fyjjz7KP//5T0aNGkVeXh7R0dG8/fbbZGZmsmTJEnbv3s3w4cMZM2ZMUPE89NBDTJ06lczMzGrXDtfu4sWLWb58OV26dGHUqFHMmzeP0aNHB/vP0Kw5ufrHUmvtUGvtYGttqrX2PqdiERGR4Hk8cMUvyvjFV9fyH65mECsAiKKEs/iYVzcez5/GzKW01OFARY6CMQ33OpyYmBgyMzPJzMw8ZEINHHLkd9SoUdx888089dRT5OTkEB4ezty5c5k4cSJhYWF07NiRk046ie+//74u/zwAh213xIgRuN1uXC4X6enpbNy4sc73ay60o6KIiByRv/4VTvjwT0zhuRqvt2Ufj68/j2fv297IkYm0fIsXL2bAgAHVzt9+++08//zzFBYWctxxx7Fq1aqgpl6Eh4fj8Xj85WCWljtcu1EVHqwICwujrKys1vZaCiXVIiIStLVr4YH7PCSxK3By0iTIyYFvviE/oSPgTazdf7me7KzQmU8p0tCWLl3K/fffzw033FDt2rp160hLS+O2227j2GOPZdWqVYwZM4YZM2ZQXl7Orl27+OqrrxgxYkSlej169CAzMxOPx0NWVhbfffed/1pERASlNfzJKZh2Q1GTWP1DRESah9/+FopKXFzNi+zsPpxbh87GvPQSuFxw/PFEv/kKnD4egJ95ZvL8H2bxy9dPczhqkeA1tefqvv76a4YOHepfUu+pp56qceWPJ554gjlz5hAWFsbAgQM588wziYyM5Ntvv2XIkCEYY3j44Yfp1KlTpSkZo0aNomfPnqSlpZGamsoxxxzjvzZlyhQGDx7MMcccwyuvvOI/f8EFF9TY7qpVob3eRJNYpzpYKSkpdvXq1U6HIU2IlsmSmqhfNIyvvoKTTvIeGwPffw/DjrHVJotuOvWX2Nmz+Ts3MzPxalZlxxMX50DAFahPSFVVl9SraUqFhI6DfaBivzDGNL91qkVEpOmrsFwtV14Jw4YBVH/6yv363xkwLJafNoRDDrz8Mlx7baOFKSLiCM2pFhGRWn3zDez/4gdclBMWBnfeeej3hrVpxXW/CYzZvPRSIwQoIuIwJdUiIlKr//51C/MYxXIG8c/jX6Z3T89h33/ZZRAW5j3+5hvIzm6EIEVEHKSkWkREDisrC3p9OJUoSujPai7LfbrWRXeTkuDkk73HKaxiztTljRCpiIhzlFSLiMhhPftkIb+0z/rLcffeVvtOFsCNQ+exkGNYxQDcz/+5IUMUEXGckmoRETmk0lLIe+5V2rEXgPwOPeCcc4Kqe+LP2nIMiwE4bs+H7M4qbKgwRUQcp6RaREQO6aOP4NIDgVHq6D/cEJgsXYu2owawIca7TFkMRSx/5qsGiVGkpZg5cybGmDqv93zVVVfx5ptvHvY9f/nLXyqVTzjhhKO61z333MOjjz56VHUPysjI4JxaflnPycnh6aef9pe3bt3KRRddVKf71jcl1SIickifPrGSkXh3WCtzRRB2zeQjqr998On+45L3P6nX2ERamunTpzN69Ghee+21Br9X1aT6m2++afB71kXVpLpLly61/uLQ2JRUi4hIjXbuhB5f/tdfLhp/LrRrd0RtxF98hv+456pPmtxudSJNRV5eHvPmzeOFF16olFQf3Izkoosuon///kyaNImDG/fdd999DB8+nNTUVKZMmULVDf0+//xzLrjgAn951qxZTJgwgdtvv53CwkLS09OZNGkSAPHx8f73Pfzww6SlpTFkyBBuv/12AJ577jmGDx/OkCFDuPDCCykoKDjs5/PGG2+QmprKkCFDGDNmDABFRUVMnjyZtLQ0hg4dypw5c6rVqzrynZqaysaNG7n99ttZt24d6enp3HLLLWzcuJHU1NTDtjtt2jQmTJjAGWecQd++fbn11ltr+SrUjZJqERGp0Zszyplk/+cvx99w5RG30f9XYygkGoA+ZavY9NWmeotPpCV55513OOOMM+jXrx9t27Zl0aJF/muLFy/miSeeYMWKFaxfv5558+YBcOONN/L999+zbNkyCgsL+eCDDyq1OW7cOFauXMmuXbsA+M9//sPkyZN56KGHiImJITMzs9L24wAff/wx77zzDgsWLGDJkiX+RHTChAl8//33LFmyhAEDBvDCCy8c9vO57777+PTTT1myZAnvvfceAP/85z8B+PHHH5k+fTpXXnklRUVFQf37PPTQQ/Tu3ZvMzEweeeSRStcO125mZiYzZszgxx9/ZMaMGWRlZQV1v6OhpFpERGr0078+pytbASiIT4IzzzziNiJaxbAi6SR/efN/Pq+3+EQazD33eFe4CeY1ZUr1+lOmVH7PPffUesvp06dz6aWXAnDppZcyffp0/7URI0bgdrtxuVykp6ezceNGAObMmcPIkSNJS0vjiy++YPnyyktXGmO4/PLLefnll8nJyeHbb7/lzFr+P549ezaTJ08mNjYWgLZt2wKwbNkyTjzxRNLS0njllVeq3auqUaNGcdVVV/Hcc89RXl4OwNy5c7n88ssB6N+/P927d2fNmjW1/tvU5nDtnnLKKbRu3Zro6GgGDhzIpk0N94u9tikXEZFqNm+GY1cEpn4waRJERBxVW7nDxsInnwJgvpkHXF33AEVakD179vDFF1+wbNkyjDGUl5djjOHhhx8GICoqyv/esLAwysrKKCoq4vrrr+eHH34gOTmZe+65p8ZR38mTJ3PuuecSHR3NxRdfTHj44VM/ay2mhiUzr7rqKt555x2GDBnCtGnTyMjIOGw7//rXv1iwYAEffvgh6enpZGZmVpueUpPw8HA8nsDmUsGMZB+u3Zr+7RqKRqpFRKSa116D+7ibB/kTO6OTif31kU/9OKjNOaP9x+5Nc+sjPJEW5c033+SKK65g06ZNbNy4kaysLHr27MncuYf+/+Vgstm+fXvy8vIO+dBely5d6NKlCw888ABXXXWV/3xERASlpaXV3n/aaafx4osv+udM793rXU4zNzeXzp07U1paWm3KSE3WrVvHyJEjue+++2jfvj1ZWVmMGTPGX3fNmjVs3ryZlJSUSvV69Ojhn/qyaNEiNmzYAEBCQgK5ubk13iuYdhuDkmoREalm+nRYQwp38iAfP70R0tOPuq1+vziWYiIB6Fmyhn2rd9ZTlCIN5J57wNrgXs8+W73+s89Wfk8t0z+mT59e6YFCgAsvvJBXX331kHUSExP51a9+RVpaGueffz7Dhw8/5HsnTZpEcnIyAwcO9J+bMmUKgwcP9j+oeNAZZ5zBeeedx7HHHkt6err/ocH777+fkSNHcuqpp9K/f//Dfj4At9xyC2lpaaSmpjJmzBiGDBnC9ddfT3l5OWlpafz85z9n2rRplUaSD37ee/fuJT09nWeeeYZ+/foB0K5dO0aNGkVqaiq33HJLpTrBtNsYTDBD8U1FSkqKXb16tdNhSBNy8KlokYrUL+pm1SoY4F1emqgo7yogrVrVrc3MhBOJyNvLXEbTd9qdjLsyue6BHgH1CamqYp9YuXIlAw52+hboxhtvZOjQoVxzzTVOh9JkHewDFfuFMWahtfbYYNvQnGoREamkwvNRnH123RNqgP9d9Tl/n+odrb7jJxhX9yZFJAjDhg0jLi6Oxx57zOlQWjwl1SIiUsl3L6+hEwlspzMTJ9ZPmyNGR8JUX/vf1U+bIlK7hQsXOh1CyNCcahER8Vu5Eqasv40tdGWu60TO6fh9vbQ7bFjgePFitAmMiLQ4SqpFRMTvg9fyOINPcGEZ5ZlLdMfW9dJur16BaSS7d0N2dr00K1JvmtMzZlK/6utrr6RaRET89r78ETF4l+ral5wGvifv68rlgkv7/MCd3M9Mzmfb0zPrpV2R+hAdHc2ePXuUWIcgay179uwhOjq6zm05NqfaGJMMvAR0AjzAs9baJ52KR0Qk1G3eDOnr3/KXoyddWK/tXxD5IWdwDwDzv+gOXHDY94s0FrfbTXZ2tn87bwkt0dHRuN3uOrfj5IOKZcAfrLWLjDEJwEJjzCxr7QoHYxIRCVnvzSjkKj70l2Muu6he2488fhjM9x63+kkPT0nTERERQc+ePZ0OQ5o5x6Z/WGu3WWsX+Y5zgZVAV6fiEREJddumfUo8+QDkdEyBChtF1Ieu5wWeVuyRsxjKy+u1fRERJzWJJfWMMT2AocCCGq5NAaYAJCUl1brXvISWvLw89QmpRv3iyO3dG8GAFYFtjneMHkHml1/W6z3Kyw2t6ERnthNrC5j97xmED+xSr/c4FPUJqUp9QmpSl37heFJtjIkH3gJ+b609UPW6tfZZ4Fnw7qioHbGkIu2SJjVRvzhyLzxdzEW87y+n3HETKUOH1vt95iek0Tl3OwBd9kQxsJG+TuoTUpX6hNSkLv3C0dU/jDEReBPqV6y1bzsZi4hIKNv0wmxa4x3X2N+2J6SnN8h99rnT/Mf5839skHuIiDjBsaTaGGOAF4CV1tq/OxWHiEioy8mBXosDq35w8UVgTIPcy5M22H8cvlJJtYi0HE6OVI8CLgfGGWMyfa+zHIxHRCQkffgh/M9O4gWuJiesLa0n1+9SehW1HhUYqW6/XUm1iLQcjs2pttbOBRpmKERERII2cyZ8wSl8wSns+vO/uH1Ew/1o6H7GAMpxEYaHroVrsfkFmLjYBrufiEhj0Y6KIiIhrKAAPv44UD7/4ogGm/oB4O4bwzpXXwBcWPZ8tbzB7iUi0pgcX/1DRESc89ln3sQaYMAA6N+/Ye9nDMzscRPZ60v4kTTuKh3IKQ17SxGRRqGkWkQkhH326m5iiKWQWCZMaJx7rh13Lc+v9x5P2IiSahFpETT9Q0QkRJWUQPq797KLJGZwCZf1/6FR7puSEjhetapRbiki0uCUVIuIhKgvZns4p+Qt4ijgEt4gpUtuo9y3YlK9enWj3FJEpMEpqRYRCVGLnp5PF7YBkBfTHjPmxEa5b8V525tX5EF5eaPcV0SkISmpFhEJQaWlkDj7TX+54LQLILxxHrPp2RP+a65iE934aXsC+Qs1B0REmj8l1SIiIWjOF5aziwO7KCb9uuE2fKkqPBx6xu6gG1kAbJ+zstHuLSLSUJRUi4iEoO+f+YHubAagIKoN5pRxjXr/fZ0G+I/zftBItYg0f0qqRURCTFkZtP70dX85b9x5EBHRqDGU9gkk1azUSLWINH9KqkVEQkzGHMu5RYGkuv31lzR6DLFDA0uAxG5b2+j3FxGpb0qqRURCzMJ/zvdP/ciPaoPrtPGNHkO7Eb39x0n71zX6/UVE6puSahGREFJWBomfzfCXc8dPgMjIRo8jeURnCogBILF8D3ZfTqPHICJSn5RUi4iEkK+/ho8Lx/IBZ1NCBB1ubPypHwCdurjYaHr5ywcWa7RaRJo3JdUiIiHkjTfgXc7nXD7gzl/uwDW+cVf9OMgY2B7fx1/ePV/zqkWkeVNSLSISIsrL4e23A+VzLm/TaBu+1ORAUmBedf4SJdUi0rwpqRYRCRFffAE7dniPO3aEUaOcjaesR2Ckujhrh4ORiIjUnZJqEZEQMfO53YRRBsDEiRAW5mw8B8ZPYAiZJHCAfw14ytlgRETqSEm1iEgIyM+HsW//lmzcPMbNXDN+k9Mh0Tm9I0sZQh4JrNNziiLSzCmpFhEJAR+8eoBzy2fSiR3czOMMcu93OiR6B6ZUs369c3GIiNQHJdUiIiFgy5NvEkMRADs6D8EMGexwRNCjh3cVEIDsbCgudjQcEZE6UVItItLCbdsGxyx/yV+OvOYKB6MJiIyEbm4PXclmjM0ge+5Gp0MSETlqSqpFRFq4D6duYCxfAlCOizY3/MLhiAIeszeRTTIZnEzJ9LecDkdE5Kg5mlQbY140xuw0xixzMg4RkZbKWih6/mV/eWvq6dCpk4MRVVbq7uk/LlultapFpPlyeqR6GnCGwzGIiLRY384t5+ydL/rL7W5qGlM/DgrvF3haMXKzkmoRab4cTaqttV8Be52MQUSkJfv2nk/pyUYA8qLaEjvxZ84GVEXckMAGMK33aF09EWm+nB6pFhGRBrJ3L/TPeMZfzrtoMsTEOBhRdR2P64kH7xIgSQWboKTE4YiaNmuhrMzpKESkJuFOB1AbY8wUYApAUlISGRkZzgYkTUpeXp76hFSjfuH1xYtwj+dDf3nj6emsamL/Lnl54bTHTTeyCMPDt9NnUNw9uQHu03z7xMaNsbzzTlcWLWqDJ2sv3diMKz6S8N7tSDuxjPHjd9K6danTYTY7zblPSMOpS78w1tr6jeZIAzCmB/CBtTa1tvempKTY1atXN3hM0nxkZGQwduxYp8OQJkb9wjuieWz/PEas+R/X8QxtBnQiecVnTodVo68jTubEsgwA9rz0Ie0uP6ve79Ec+0RBAfzhD/CvfwXOPc11XIf3hAfDj6QxM/xiSidN5qZHu9K+vUPBNkPNsU9Iw6vYL4wxC621xwZbV9M/RERaoK+/hkVr4vkX1zEqbgmJn73hdEiHtCcx8LDivkUbHIyk6dixA447rnJCDbCdwMotLixDWMo9ZXdx1397M8P9B174225NDxFxiNNL6k0HvgVSjDHZxphrnIxHRKSleOyxwPEvJhkS3K2dC6YWJZ26+4+LVm9yMJKmYfduOP3kEuyPP/rPnXEGzJoFd0zrhx0+nJI+A/C4wvzXoynmhuK/c/HtvZja5wk2rPM4EbpISHN69Y+J1trO1toIa63bWvuCk/GIiLQEK1fCe+8Fyr//vXOxBMN2CyTVJiu0k+rycrj8F+XcvvIK5nMcp5lZ/Pvf8NFHMH48RF45EfPdd0T+tAJXzj7sK6+yr+8If/1W5PL7TTexNWUs7z650blPRCQEafqHiEgLM/O2+Yzjc8By3nkwYIDTER1eRL8eHCCBH0klm/p/SLE5+etfYeSs+7mUGcRRwEeus5kyajnG1PDmhATMLybSZvV8SmbMZFdS4AudVp7Jb37v4vrrtaCKSGNRUi0i0oJs3QonfnArnzOe7xjBvecvdjqkWkWOO5HW7GcwP/JY50edDscxq1fD5/d8zV3c7z8Xdt21MHDg4SsaQ+Ql55OUncmWa+6mjDBu4nGy6MYzz8C4cbB9ewMHLyJKqkVEWpJ3b53HifZrAIaaTNLHN/3lILr3MOBbq3pTiM7+sBZuur6Yf5dfQxje+dB2zEnwxBPUPExdg8hIuj5/L8U/LCP/kqv9p+fNg2HD4LuvihoidBHxUVItItJC7N0L3Wf8zV/eevJlkNz0p1N06xY43rzZm2CGms8+g9QvnqQfPwFQHt8a88rLEBZWS83q4ob1Z/prhocfBpfvp3zbrT+SfFJPZl//dn2GLSIVKKkWEWkhXvndd5xV9r6/7H7qVgejCV7r1t4XQFER7NrlbDyNzVp44o5dlad9PHgfuN1H3aYxcMst8Mkn0L3NAd7kIjqznfHPXMis9D9Skq/NYkTqm5JqEZEWIGuzZegrf/SXs4+7CNegJv6EYgWjO/7Ez3mNW/kbe9+f53Q4jWrWLBi98AkSyAOgtO8AuO66emn71FNh7hvbiI8IPK146pLHWNl5HFt/2Fov9xARLyXVIiItwMyr32e0by51mQmn63//4nBER+Yy+xKvMZG/cTv240+cDqdRPf/IPn7DP/zliAfvhYiIemvffUoKiesWsrDLOf5zQ3LnEj5iKD88Mqfe7iMS6pRUi4g0c8syyzj189v85e3nX4fp19fBiI5cedfAWtV240bnAmlka9dCz9nP0opcAEp694cJE+r9PnHJbTlm87t8feZfKPf96O9gdzL01vF8cfL9lBdrG0aRulJSLSLSjFkLn098ngGsAiA/vBXuf9/lcFRHLrxPD/9x1LbQWQLk6afBg4u9tAEg8o5bj+rhxGCYMBcnfvR/LHt8NrtcHQAIw8O4jLtZ1XEMW79e1yD3FQkVSqpFRJqxN5/ayhWr/s9f3v/r2yEpycGIjk78oMBIdcL9hgoVAAAgAElEQVS+0EiqS0rgpZfgUW6hK1tY+of/wqWXNvh9h/z+ZFi0mKWtT/SfG7T/W2LGDOfVZ/aH5OorIvVBSbWISDO1cyfcdG8iL3ANHgy7W/eiy8NNfE/yQ2iXHlj6r11hNpS2/NUpPv0U9uzxHrd3x5D68BUQE9Mo904a0oVBO+cwa+yDlBIOwCP8kUnXt+b887VZjMjRUFItItJM/fa3sGVfLLfwKBd2+obYt15utKSsvnXrF802OgHeKQls2eJwRA3vf/8LHE+aFFhTurGERYZx6pw/sez5BbyRMJm/4Z2X/9573k0cn34ayks9jRuUSDOmpFpEpBl66y2YMSNQ/vW044g95XjnAqqjTp1gswlMASlc1bKngBw4AO+9G5hncfnlzsUy9JpjOGvbi1x3Q2Au9759MPWGFWxKGMSyB2aG5o48IkdISbWISDOzYW05U64OrNZw+eVw+ukOBlQPXC7YFdvDX963eKNjsTSGDz+E90pOZzqX8tse7zGob0ntlRpQXBxMnQqzZ0Pv3gCWp7meXsWrSL1rAj8lHsu6f3yk5FrkMJRUi4g0I/n5kDHqT7x54FQ6sp3u3eHJJ52Oqn7ktguMVOevaNkj1V+9spnTmMWlzODxzRNg/36nQwLglFNg2TKY+vt1DGap/3zfA4vo/duzWdn2BJbc/x62rNzBKEWaJiXVIiLNRHk5/GfMf5i882FOJoPFDOWdv6+nTRunI6sfpZ0DSXX5hpabVBcVQdysd/3l/OPHN6kVW6Kj4YbH+1CweA0fD/ojBQTm6Q/Imc+Qu3/Glri+LJz4KIVb9joYqUjToqRaRKQZsBaePftdfr3oV/5zpWnDSP9Z98PUal7K+qfyCafzb6awpP14p8NpMJ9/DmeWvOMvx08638FoDs2d3p4zlz1C1px1fNj7NxQTGbhWsoFhr92CdbuZM+hG5s3TzBARJdUiIk1cWRk8c+rbXP3pJYTj/bP71vaD6TZveoNtFOIEe+IYzuQTfs2/+SBhotPhNJgv3trHSXzpL5ufnedgNLVLGduZs9c+xdYv1/LZ0FvZQ1v/tVgK+WlFCaNHQ58+cPPNMGcOlJYow5bQo6RaRKQJKyyE/xzzD379+UVE4X2YbUd8bzplfgIJCQ5HV7+6Vxh039RCZ39YC+Xvf+T/5ehA/xHQpYvDUQWn55hkTlv0N+zmbD644AVWRKYDMB3vL0Dr18Pjj8O4cfB53Hks7XQq31z0dza8tQhbqm3QpeULdzoAERGp2YZVxXx38m38anvgScTtCX1ot2g2rq6dHYysYXTrFjjevNm5OBrSTz/B8N0f+cuxE5v2KHVN2ifHcM7bV2M9k8l8aQn95qax6A3vMoEAseRzctlnRO0ogbdmw1uQb+JY1/449qeOInbMMLqfl077oclgjLOfjEg9UlItItIEffCXpXS/63J+7gmswLCp80i6Zb6P6dB0HmqrTxWT6uxs77SX8Bb2U+rTjz38nFn+cvg5ZzoYTd0YlyH9qnT+fRU8ORUyMuD992H3G98RtavyEoFxNp/Buz6HOZ/DHOBeyDGJZLUZwmtXfkzv1BhSUiAlBdq3d+KzEam7FvbtSkSkeVu6FG67DQZ8Mpu/V1jSbG3qz+iz4FWIjXUwuoYVHQ2TW73F0AMZuMuz2TfzdyRdPNbpsOrV2jcz6cAuAArik4hNT3c4ovoRHQ1nnOF92akns2JWFlte+IToebPptW0eXT3Z1eok2hwK967mL49X3gX07ISveLzsN+xp24eSDsmYbslE9U2m9aBkOgxLJnFAZ0x4y3mWQFoOJdUiIk3AwoXwxBPwyiveebef8nvO4z2OMwvYftPD9Hnkhsbfx9oBZ0XO4iL+DcDG+eNaVFJdVgaJ333mL5effGqL/JoaAwNPczPwtF8Cv8TjgVWzN7PtzXkwfz6tNy6hZ+4S2pDDEoZUq5+cu5y+LKXvlqWwBVhc+XoZYex2dSQnqgNLOp/JJ2P+QseO+F9ditbT3u4itksiCcmJJPZsQ0RcZLX7iNQ3JdUiIg7JXpXH0sc/J+qNl3l731he5gb/NeNy8dml/2X4zYX0GNbfwSgbV1H7ZNjtPS5el+VsMPVs0SLoUrIBDwYXlvgLm/k2mEFyuaD/ad3of1o38D3UWFZq+enLzZQtLuDuPFi92vtaswb6Fvx02PbCKaeTZyudCrfy/fpUpq2vfP1hnuEWHq10roAYcl2J5EUkUhDVhuLoRDyuKJ4fspufhlxEfDzEx3t3luy5ZS6ty/YQkRhHVNs4ItvGE9M+jpikeGKT4ohqE4txaS64VOdoUm2MOQN4EggDnrfWPuRkPCIiDaW0xJK1aBeb3llMYcYCElfMY1huBmf5VvTozHKe5nrAcPbZ8NBDkJractagDpanixtW+QpZ1acMNGdffgm38m9u5yHuGjuXm88+3umQHBMeYeg7vjt9x8M5Fc57PLB1+V0smXcxuT9upGR9FiYri+hdWbQ6kEWH4iyS7C7/+3fSoVrbHdhZ7VwshcR6CulYvA2KAd9DlV9tH8zDn15U6b0f8FdG8FG1NvwxYigkmmITzZ2t/8HHbSYRHY3/9aeNU+hUmoUnIoryiGg8kd6XjYqCqGhsVDSeqGiIjGLDMRdS1LE7EREQGQkREdBtyftEhHlwRUXgioogLDrwCo8JHEfERmCS3UTERRIR4V1dM8xlCTflhEWG4Qozeg60kTmWVBtjwoB/AqcC2cD3xpj3rLUrnIpJRKQ21kJJkYei/cUU7y+iZH8hJQe8H7fH9WbH/mh274bdu6Fw8y7O+vhG2u5dS3LxWnpxgF6HaHcgK7njrEzOvnMox4durkV4z2T/ccSOljVSnZHh/ZhDGxIvPxf0QF41Lhe409rgTjseqPl/hIK9RexZuZP9a3YwpKQtz7pgxw7va+dOiPg+mRU7jiW2NIf4shxa2xwiqHlJvxwSq52LI//wMWK9SbotZH+OZUNO5eu9mMtAVgb1+T7wwRC+oPIvz/u4nESC27Z+EMtYwSB/OYZCCogDoBwX5YRVenlM5Y+ndF5JYUQrwsK8DwW7yzfxTNY5eFze93hMGNZ3bE0YHpev7AqnJDyWR0a+hcuF/+XOW8X5Kx70zgEyLqzLhTEG6zvGuLzXXC5y4zoya+RdB4u4XNB1VyZDV79GpZO+Y+PytnPwOLdtd1YNm1TxLXTcsoiuG+YG6oW5MP4AvTEdPM7v2Ivd/U/0huqr32ZTJvu/X8hXn5eyN3VMUF+DipwcqR4BrLXWrgcwxrwG/Aw4ZFIduWYTGyJTgmp8QveF5Jt4fzm5dD0vbAnuT22FJo7zumVWOjekaAGP7ZhU+Y2HWNs+K6InV3WZVencKfnvccfu3wd1/yVRI/l9x+mVdqe6JPd5bsh5MKj6s2PO4752T1Y6d2POA1ya+1ylc4damv+1+F8xNfHOSufu3Xsj4wveC+r+/2h9J9Pjp1SK/5ndF3NMybdB1b8rcSqfxlTeYeytnSfSrWx9tfemANurnLu27RssjDqh0rmvt/cl1uYFdf8L2s9lU3hvf/yxnjy+2dUnqLoAx7dfS4Er0Pd6lK3l/b0nVH9jDV+AAhPHiKQNlc4NL5nHtJzzD/n1qmhTWG/Oajvf27yvwjnFb/Jo7rVBxf5dxCgub1X563x14T+5veCuoOp/EHkhN8VX7mf/V3AnvyyaGlT9F6Ju4MFYbz8/GP9T+ddwXsmbQdV/IOYBnov6DWVlo/yrRryedxbHlX0dVP0bY1/k3YiLK537Ojednp61ABgs4ZQRRQlRNdS/lMUsIfDgWQxxPMTrh73n+vg09h1/Nr3uvowHRg867HtDQWw/t/84PqfljFSXlcHXFbrhSSc5F0tzF9s2mthR3Uge1Y3UGt/xgO/lZT2WAzsKOLBpH3nZORRuy6Foew6blqxgePrp/DUe8vMhL8/72vfNaBbsaUVEcT6RpXlEleYTXZ5HjCefOJtHDEX+touIrnb36ArXa1NKRLVzEZQedf2KdcPwEIYHKrZnK3/M2mKo/JOxkL4sC+reucTzzjuVz41lG4/xclD1V9OPcxdU/tkykRXcxN+Cqj+HsdzyWuW87A98waXcElT9V5nIZZxY6dxjvMTNPA7AenoG1U5FTibVXYGKwxDZwMjDVYiimJ6la4JqfO1aW6mjGErpwdqg6uYSz9oqb+1KId1ZF1T9orIw1lV56wjy6MaGmitUsb6sG+ur5I+l7CeZjUHVj8rdxcbcqmf30ZXgFn41+/exqcovyZHsoSvBjRqV7s1l897K5+LYRWe2BFU/f08hVX+UJrKDTmwNqv7+3SXV7tSe7SQQXFK9e2d5pTvFY+nAjqDqAuzcWbnvtaKc9uw65PsryrWF7Khyq1xKaXdwkmkt9pS1ZWeVv3wWUUJb9tZcoYrokgPsrnKrMopow76g6ocV57OnuOrZQloHOepii4rYW+XnURgFtDr4t9palBeWkFMIVPhBE0FB0F/7koKyapFGUUB8LSNXB7VjT6VyIbFk0xW3r0fmmXiyEwawv98Iok8aSa/JJ9FrULeamgpZbdICSXWbgi3e+QAt4GG+JUsg1/d9uWtX6HWoP1lIvTMuQ6vOcbTqHAcE+ldxRjxjxw6oocYDNZwLKC8uo3h/EUX7i3nSFcffgKKiwGvvD/9jd84ByvOLvK+CImxBEZ6iYigswhQXYUqKcJUUMy41mUGxUFoKJSXej4vnnkdEaQGu8lLCyku8Hz2luDzej+GeUsKs97hV22jalnt/aSsvh/iyMsqLXb5kunblVF5FJcy3KdHR1D3S+p4a9h90BRl3fdS3VJ8bYyoMX33OKcDzQbcHzibVNc30qTYYZ4yZAkwBGNbQEYmIBOngnErvK4oSVxR9knMpc++idetS36uEOXvvJLFbBLGD2xPWuZV/s4tCYOGu9ZBR/S8woSx7Twx7aEs79hJhS/nmnXcoadu29oq1yMvLI+Pg/AsHvPdGEnfxFl8xhug+Pfnyy+AGWaThNFSfKB0CEON7BVRNAT3AWDYBlbcPLbv22kNMVqnuUdYDlb+HfM3nvq07PXhKPXhKLZ5SD7bMe2zLAtdebJ2Jx7rweAzl5QZTXMLbW1/BlllsmQfKPf7jg2XKyqHcUo6Le3svw1qw1uDxQNyBCN5Y/1fwWLAWYy14PN6/PHo8GN95PJaCyFbc0O8nrDX+Njrtbsfbm27x/jJtwVgPxnqwHu8xFrDedrbHdeOSXlnet/raiN2RzAdbJ/vqVYjBenD5yx4MHna3Hsy4Lt5RrIP1i7f1JGPXWaTk/8jSfufAj0eWVBtrg/mjcv0zxhwP3GOtPd1X/j8Aa+1fD1Wnj7u3/ey5Qz88UFFpj77eeTcHU/eSEiK21r7vrTHe/5R2r/znflNYQPjO6iOlNT0EYMMjKOtaeY6UycslbO+uoB4asFHRlHfsUql914EcXPtrHy00BjwxcXjaV354w5WzF1d+teHrmu8fn4AnsfIPMdfe3ZjCgqDi97RKxCa0qhz/7p2Y0pLD1KoQf2JbbEzltXhdO7djPNV/A164aCHDjqn865anbXuIqvzHedeObZXiOWz87ZL8O04YA3g8uHZXf/DlkPG371B5ZK2sDFdOzSPFVeOxGGz7Kht7lJRgDuwPLv6wMGybwNfOGKC4GJMf3EitiQjHtmpd+WRhIaaosOb3V40/ItL7CH1FBQX+r31t8dvIKIgJ/CAyBsjPx5QH+SMmOhqiopg7dy6jR4/2nsvP936DDiJ+oqO9TwpVlJ/vn4sSFQVRceG4YqK0E1wDKSiANXHppLMEAM/873CNHF7ndjMyMhg7dmyd2zlad5/6DffNHgXAnk4DabdtuWOxiJfTfUKapoyMDMb65mcZl2uhtfbYYOs6OVL9PdDXGNMT70qUlwK/OFyFsLhwep0Z3Jzq6iIhte9R1gWIBYKfV1tdgu91tBJ9r6PV1vc6Sj3q+ERNt+pPaB+R5E41nl6T/xNdhnetvb67Lls6u8Bd8/2DEw7Jdfn8I4G67KAX5XsdreojLkcm1vc6Sm3jjrhKQkIZbdr4Cm2OvH4liXWsL0ckNhZ2RLih1JtUH1iRTWI9JNVOi/hhfqAwLOif0SLihKMcNHFsopq1tgy4EfgUWAm8bq3Vr+4iIiFuf0JgBZADK5r/CiDZ2dAvZ4G/nHj6YR8fEpFmytF1qq21H8FhFoMUEZGQs6zXeSzd25Vs3FzW8wSa+6Oc334LIwkk1WGjjnMwGhFpKNpRUUREmpTtQ8/kuR/OBGC4hfEOx1NXyz/fzsW+h9FKw6OJSEtzOCIRaQi1Tv8wxsQaY+4yxjznK/c1xpxTWz0REZGj4Q6sekZ2C1iqunRuYJT6QJ9h1R+GFZEWIZg51f/Bu6nnwa2NsqltEUcREZGj1JKSao8H2vwUSKojT9R8apGWKpikure19mF8W/JYawupeY1pERGROmtJSfXatTC0JJBUx4/XfGqRliqYpLrEGBODb2MWY0xvvCPXIiIi9c7thheZzPccy+tfd4Jt25wO6agt+r6c4XzvL5vjNFIt0lIF86Din4FPgGRjzCvAKOCqhgxKRERCl9sNxWQylEwoB7s5C9O5LmvNOyd71kpa4d14KzeuEwnJybXUEJHmqtak2lo7yxizCDgO77SP31lrdzd4ZCIiEpJatYLtYW4ozwQgb1U2CSNHOBzV0Zm/pi238jdGsoBjTmhPgnbiFGmxDplUG2OOqXLq4N/fuhljullrFzVcWCIiEspyEpIhx3t8YHlWnfajdYrHA7OWd+EAtwKw6XmHAxKRBnW4kerHfB+jgWOBJXhHqgcDC4DRDRuaiIiEqsL2bn9SXbSueT6tuG4dHDjgPW7fHjTzQ6RlO+SDitbak621JwObgGOstcdaa4cBQ4G1jRWgiIiEnvLOgQzUbm6eW5UvXBg4HjYMNPNDpGULZvWP/tbaHw8WrLXLgPSGC0lEREJdeI/AunoR25vnSHXVpFpEWrZgVv9YaYx5HngZ77J6lwErGzQqEREJaTH9AiPVcfua50h12w//x4/8jUzS6e25BDjP6ZBEpAEFk1RPBq4DfucrfwU802ARiYhIyEsc1DVwXLAVysshLMzBiI6MtdBm3Q+kspxUlpNTkoKSapGWLZgl9YqAx30vERGRBteldwy7aE8Suwm3ZbBjB3Tp4nRYQVu/HgaUZPrLrU/SrEmRlq7WpNoYswHfbooVWWt7NUhEIiIS8txu2EgySXi3RbCbszDNKKlekmkZxxJ/2QxVUi3S0gUz/ePYCsfRwMVA24YJR0REBNq0gcsiH6G4BLJx813vnrR2OqgjkDV3E4nsByA/ui1xbnctNUSkuQtm+seeKqeeMMbMBe5umJBERCTUGQNru5/CTz95y1t2QeskZ2M6EsULAlM/9vcYQpzW0xNp8YKZ/lFxZ0UX3pHr5ri5lYiINCNuN/6kOjsbBg50Np4jEbsmkFS7NPVDJCQEM/3jsQrHZcAG4JKGCUdERMSr4oyJ7Ga0VHVREbj3BOZTtzlZSbVIKAgmqb7GWru+4gljTM8GikdERAQIJNWRFJOzci/Q2dF4grVyJQwhMFIdNVJJtUgoCGZHxTeDPCciIlJvUuKy2U5HionmqmdGOh1O0FYvyKEnGwEoNRHQv7+zAYlIozjkSLUxpj8wCGhtjJlQ4VIrvKuAiIiINJh2Ke3pyE4AWjejDWD2frXMf7y7wyA6R0Y6GI2INJbDTf9IAc4BEoFzK5zPBX7VkEGJiIh06RXNTpLowC7CbHmz2QDm/X2juY9tDGEJd15V1kwmrYhIXR0yqbbWvgu8a4w53lr7bX3e1BhzMXAPMAAYYa39oT7bFxGR5s/t9q5R3YFd3hNZWc0iqV62DHbQic/oxONXOB2NiDSWQ86pNsbc6jv8hTHmqaqvOt53GTAB+KqO7YiISAvVvj1sNYElQIrWNv0lQPbtC6xUEhkJffs6G4+INJ7DTf9Y6ftY76PI1tqVAEaL4YuIyCG4XLAvPtk76RA4sCK7yT/QsywwnZoBAyAiwrlYRKRxHW76x/u+j/9tvHBEREQCCtq5/Ul1YTMYqV79Qy7HsIaVDCA1NdbpcESkER1u9Y/3AXuo69ba8w7XsDFmNtCphkt3+OZrB8UYMwWYApCUlERGRkawVSUE5OXlqU9INeoXLcfemMAGvjk/rmDDUX5dG6tPZL2+lYVMwoNh8Tenk5FxW4PfU46Ovk9ITerSLw43/ePRo2rRx1o7vi71K7TzLPAsQEpKih07dmx9NCstREZGBuoTUpX6RcuxerD1T0ZsW7CfIUf5dW2sPvHjzscBcGFJ6pfMMPXDJkvfJ6QmdekXh5v+8eXBY2NMJNAf78j1amttyVHdTURE5AhE9wk8qBi7t+lP/2izZbn/OH7kIAcjEZHGVuuOisaYs4F1wFPAVGCtMebMutzUGHOBMSYbOB740BjzaV3aExGRlqnVwEBS7SopAnvIWYmO27MHehcHkurEUUqqRULJ4aZ/HPQYcLK1di2AMaY38CHw8dHe1Fo7E5h5tPVFRCQ0dOkdQ1/WsJUu9B8Ux8ImvGjU6lWWQazwl12pAx2MRkQaW60j1cDOgwm1z3rw7RsrIiLSgNxuWEtfCojzr//cVGXN30JrDgCQF5EInbWXokgoCWakerkx5iPgdbxzqi8GvjfGTACw1r7dgPGJiEgI69QJwsKgvBx27oTiYoiKcjqqmuUtCEz92NNxEPHai0EkpAQzUh0N7ABOAsYCu4C2wLnAOQ0WmYiIhLywsMoDvlu3OhdLbVwrA0l1SV/NpxYJNbWOVFtrJzdGICIiIjXp0aWEyOws3GSz74tW9LxmqNMh1ahVdiCpjhqq+dQioabWpNoY0xP4DdCj4vtr2/xFRESkPvzcM50buQqAjS9OhGtedTagGpSXg3t/IKluN0Yj1SKhJpg51e8ALwDvA56GDUdERKQy0y0ZfvAeh21rmk8rbtwIWdZNV7Jxs4W4EUqqRUJNMEl1kbX2qQaPREREpAYVN4CJ2dM0k+rVq+Fi3gTgrBNy+LBTa4cjEpHGFkxS/aQx5s/AZ0DxwZPW2kUNFpWIiIhPqwFd/cetc7PB4wFXMM/ZN57VqwPHXQclghb+EAk5wSTVacDlwDgC0z+srywiItKgOveJYy9taMs+Imwp7NoFHTs6HVYlFZPqlBTn4hAR5wSTVF8A9LLWljR0MCIiIlW53ZCNm7bs857IzlZSLSJNTjB/P1sCJDZ0ICIiIjXp0sWbVB9UvjHLwWhqdvzip/k1zzCWOfTvVuB0OCLigGBGqjsCq4wx3xOYU22ttT9ruLBERES8IiNhd0wyFHrLB1Zm08bZkCrJzYXr9v+VZLwPUZaFrQI0XC0SaoJJqv9c4dgAo4GJDROOiIhIdflt3P6kumBN00qq1y7JZ6gvoS4lnIh+vRyOSEScUOv0D2vtl8B+4GxgGnAK8K+GDUtERCSgtENg+kfZxqa1rN72r9b4j3fG9YKICAejERGnHHKk2hjTD7gU76j0HmAGYKy1JzdSbCIiIgCYZDd7MtuSRTKlEW66Ox1QBXkLA08p5nTuT9fDvFdEWq7DTf9YBXwNnGutXQtgjLmpUaISERGpIO+48bR/fw8AfxgKwx2OpyKzepX/uLy35lKLhKrDTf+4ENgOzDHGPGeMOQUtZy8iIg5wJwd+/GzZ4mAgNYjfEhipjklXUi0Sqg6ZVFtrZ1prfw70BzKAm4COxphnjDGnNVJ8IiIiuANTqsluQlOqrYVO+wNJdftRSqpFQlUwDyrmW2tfsdaeA7iBTOD2Bo9MRETEp6km1VuyLX1s4EHFxJFKqkVCVTBL6vlZa/cC//a9REREGkXXrpDGUgazlO6bs/FknoUrfbDTYbFx3hZGkw/A/vC2tE5q73BEIuKUI0qqRUREnBAbC3+M/AdXlDwPHjjwaWtaNYGkes83gakfO9uk0Nro0SORUKWkWkREmoW8RDfs9B7nr86ilbPhALAsx82P3EF/VtEpbRB9nQ5IRByjpFpERJqFkg6BpLpsQ9OYWP31zhQ+5QEA3rze4WBExFG1PqgoIiLSJCQn+w/N1qaRVK8JPKNIip5RFAlpjiTVxphHjDGrjDFLjTEzjTGJTsQhIiLNR2SvwBIg0buyHIzEq6gINm70HhsDffo4Go6IOMypkepZQKq1djCwBvg/h+IQEZFmImFAIKludSDbu0i0g9auDYTQowdERzsajog4zJGk2lr7mbW2zFecj3f9axERkUPq2LcV+32PJ0aWF8HevY7Gs23WMuYzkmlcye9jn3U0FhFxXlOYU3018LHTQYiISNPmdkM2TWcXmLwFyxnJd1zJS5yS/56jsYiI8xps9Q9jzGygUw2X7rDWvut7zx1AGfDKYdqZAkwBSEpKIiMjo/6DlWYrLy9PfUKqUb9omfLywojAzSBWALD0w4/Yu29fkHXrv0/k/rDQf7yzbUd2qc81K/o+ITWpS79osKTaWjv+cNeNMVcC5wCnWHvoiXHW2meBZwFSUlLs2LFj6zNMaeYyMjJQn5Cq1C9aJmvhf+Fu71AM0DuqDYOD/Do3RJ/4NOc5/3Gnk45jgPpcs6LvE1KTuvQLR9apNsacAdwGnGStLXAiBhERaV6MgXXtRvD+jp1k4+a0xP70digWa6FDTmA3xXYnaD09kVDn1OYvU4EoYJbxbuk631r7a4diERGRZuKbtGu5b8e1AHzYGceS6t27LH3KA0l10mgl1SKhzpGk2lqr1TxFROSIuSs8p7hli3NxbJi3lRHkAXAgLJFWHTs4F4yINAlNYfUPERGRoLibyOIfu+cFRql3JKZ456aISEhTUi0iIs1GU0mqi5YEkur/b+/eo62s7zuPv7+A3PFgDihXuQkH8FKIxIipBipxRaukNk5rmsbm0jF24ujM2JmkcWZik7a56NJkdaVt0sEdhpAAABRQSURBVKQ2Nkmtk2RSL4mTgD1JUYwY5SaCXEQEvCAieOTmOfzmj7M9e3s8wMZ9O8/D+7XWWT6/vZ/9e74sv+s5Hx5++3leG+fSD0mNW1MtSdIxGzsWruerTGU977xnK+z6RzjppLrX0W9TMVTH9Ol1P76k3sdQLUnKjHHj4JN8gxmshR3As882JFQPf6EYqofN8Uq1JJd/SJIypDc8VbG9HT6x/2+4lHu4gVs4+dJz6l6DpN7HUC1JyozmZtjeZ3zXeN+6LXWv4emnYUPHJO7jUu4ccwNDp487+ock5Z6hWpKUGRGwu+nUrvFra56pew3riis/aHHlh6QCQ7UkKVP2njyha7t9o6FaUu9gqJYkZUrHuIld2/Fs/UP1syt3cQIHAUO1pCJDtSQpU06YWrxSPeiF+ofqS+6/jr0MZh3TeM+en9b9+JJ6J0O1JClThp8+jkN0PsFw6Kvb4eDBuh5/5Mvr6EcH01jP2NMG1fXYknovQ7UkKVPGT+nPdsYA0IdU19vq7X4lMbm9uKj6lAtc/yGpk6FakpQpEybAMxSXgPBM/ZaAbFr6Ak3sAeDVPifSd+youh1bUu/mExUlSZkyYQLcysf4CZewte9Ebm+ZUbcrRC8tWdu1/XxTC8Mi6nRkSb2doVqSlClDhsCPm/+YnTuBDvgiFBaD1N7ex4tLP9rGuPRDUpHLPyRJmTOhMas/6LuhGKrTNEO1pCJDtSQpcxoVqk98vhiqh5xtqJZUZKiWJGXOqcUnldctVB86BOPaimuqR3nnD0klDNWSpMyZOL6DH3E5jzGbP/nzUzoTb41t3bCfCWkzAIcImt41rebHlJQdhmpJUuacOqkvv8kSZrOcE/e9CM89V/Njbnl4Oy/zDgCeGzAJBg6s+TElZYehWpKUOY24V/XyPZM5mR2MYAd/e8k9NT+epGwxVEuSMqd7qE6bax+q1xW+o7iTEQx798yaH09SthiqJUmZ09wM2/pN7BrvX7u55sdcV7zxBy1+R1FSN4ZqSVLmREBbc/FK9Wtr6nelGgzVkt6qIaE6Ir4QESsjYnlE/Cwi6vUwLElSTrw+phiqO56ubaje99ohfmvL7ZzLUpr77GLKlJoeTlIGNepK9c0ppbNSSrOAe4H/3aA6JEkZ1WdSMVSfsK22oXrzkq3czsdZynms5zT696/p4SRlUENCdUppT8lwCJAaUYckKbuGzCyG6iE7n4FUu18lL/yi+NCX7U0zanYcSdnVr1EHjoi/BK4CdgPzj7Df1cDVACNHjqS1tbUu9Skb2tra7Am9hX1xfNjZMZLdnEgTexjQvpcHf/xjXj/ppB73rbQnti9e0rX9wvBx7LC/Ms/zhHpSSV9EqtHf7CNiETCqh7duTCn9a8l+fwYMTCl97mhztrS0pHWl3xTRca+1tZV58+Y1ugz1MvbF8WHJEhhy/mxms7zzhaVL4dxze9y30p74+Wl/wvs2/h0Av/oPt/Duu25423Opd/A8oZ6U9kVE/DqlNKfcz9bsSnVKaUGZu34fuA84aqiWJOkNEybAx/kKHfTllXdM4bF3javZsYY/X1z+ceI502t2HEnZ1ZDlHxExNaW0vjBcCKw90v6SJHU3Zgy09nsf7e3Ay7D3AAweXP3jdHTAuNeKv6ZGzzdUS3qrRt3940sRsToiVgIXAdc3qA5JUkb17Qunnlocb95cm+NsXf0Ko3kegP0MYPisibU5kKRMa9TdPz6YUjqjcFu9y1JK2xpRhyQp2yZPLm5v3FibY2xbXLxKvW3ItM40L0nd+ERFSVJmvfEQlkHsZdfS2qwk3PNIcd5XTnHph6SeGaolSZk1ZWIHWxjPXoZw1RdnwL591T/Ik092bbZP9R7VknpmqJYkZdak0/rSXvqd+xosrF627wzu5jI2MYmB75xZ9fkl5UPDHv4iSVKlJk+GTUxmEps7X9i0CWZU92ry117+CDv5CADPXFPVqSXliFeqJUmZNWUKbGRK1/jQ+g1VnX/HDti5s3N78GAYV7tbYUvKOEO1JCmzmppg26CpXeO9y9cfYe9jV7KcmunToY+/NSUdhqcHSVKmtY2Z1rX9+hNPVXXutSU3FJnujT8kHYFrqiVJ2TZ1KhTuUd1vc3WvVI/83lf5Jk+wmjOYOmohMKmq80vKD0O1JCnThv7GFA7dH/QhMeSlZ2D/fhg4sCpzT1p9N5fzbwA8dMIkDNWSDsflH5KkTJs8YwCbmQhAH1LnHUCqISXGv7Kqazhi/pnVmVdSLhmqJUmZNm0arKf4ZUWeqs666j0bXqT50EsAtDGEie+dUJV5JeWToVqSlGlTp8JTdH5ZcVuMJe2tzlMVt/xkddf204NOp/9Af2VKOjzPEJKkTGtuhq823cQQ2hiXtrL9vR+qyryvLCmG6p1jzqjKnJLyy1AtScq0CBjR0sxehgCwvko3AInVxfXUh2YYqiUdmaFakpR504q3qq7WkmqGby1eqR52nl9SlHRkhmpJUuZNLfmeYjWuVKeOQ5za9kTXePzFXqmWdGSGaklS5k2bBiN5kfk8wJT7vw7r1lU037alWxhGGwA7o5lTzjqlGmVKyjFDtSQp86ZPh5v57zzAhVyz+lp44IGK5tv+0xVd2882nUH0iUpLlJRzPlFRkpR5LS1wJzO7xu0r11T0C+6X6Xy+wN3M5nGmzxnPrMpLlJRzhmpJUuYNGgQvj5oJz3eO9z22hmEVzPfIhndwL5dxL5fxrSurUqKknHP5hyQpH2YWr1T3e2pNRVOtXFncPuusiqaSdJwwVEuScqH57InsYyAAg155HnbufFvz7NtXvINIBJx+erUqlJRnhmpJUi7MOKMvT1CSgFesOPzOR7B22av0O3QA6LxV3+DB1ahOUt41NFRHxJ9GRIqIEY2sQ5KUfTNnwvLSrxQuX/625jn0tb+mjaE8xmw+1fTdKlUnKe8aFqojYjzwPmBLo2qQJOXHjBlvDtUdj729UN1/+SOcQDuzWc7E0QeqVZ6knGvklerbgP8BpAbWIEnKiSFD4MXRxVB98JG3sfwjJUY/+6uu4YkXnVuN0iQdBxoSqiNiIbAtpfT2FrxJktSD/nOKt+oYsGkNHDi2K83tm7Yw4vXO+/Lt5kRm/O6MqtYnKb9qdp/qiFgEjOrhrRuBzwIXlTnP1cDVACNHjqS1tbVaJSoH2tra7Am9hX1x/Oo/YgKLuJBdnMTrLZMYv3gxHYMHl90THf/8EBcWtleccDaH1v2SJyt74rl6Kc8T6kklfVGzUJ1SWtDT6xFxJjAJWBERAOOAxyLinJTS8z3M803gmwAtLS1p3rx5tSpZGdTa2oo9oe7si+PXq6/C+25fBMAFI+AXl3S+Xm5PrPrLu7u2X5x8HlfYR7nleUI9qaQv6v5ExZTSKuDkN8YRsRmYk1J6qd61SJLyZfbs4vby5ZBS572myzV41cPFwbvfXb3CJOWe96mWJOXG2LHQ3Ny5vWcPPP30MXx4717Gv/jrruHJlxmqJZWv4aE6pTTRq9SSpGqIgFklt6p+/PHyP7t38VL6p4MArGEGZy04+SifkKSihodqSZKq6eyz4Xf5Id/hKi74j9Pg5z8v63PP3dnatb2yeT7Dh9eoQEm5ZKiWJOXK3LlwIYu5in9i5K71sHRpWZ/bvK0/2xkNwGtz5tWwQkl5ZKiWJOXK3LmwlLld40O/XFLW5/4i/hdj2cY01jH8yvfXqjxJOWWoliTlyimnwObxFxRfWPLv9Dl48IifOXgQHn4YIFjPNOZeNKymNUrKH0O1JCl3Jr53AmtpAaDPgf00rVp1xP2XLYP9+zu3J0+GMWNqXaGkvDFUS5Jy57zz4GclD+49admyI+5///3F7fnza1WVpDwzVEuScuf887uH6kcPv/OKFVz015dyNd9gDNu4+OI6FCgpdwzVkqTcOf10WDdqHgc5AYBhmzYe9kkwbXf8iPN338c3uIab49MsWFDPSiXlhaFakpQ7EXDBJUNZRElC/t733rpjRwd85x+7hmun/w5NTbWvT1L+GKolSbl08cVwB1cVX7jjDkjpzTstWsTQnVsAeIlmhn7osjpWKClPDNWSpFxasADu7fMB9jCMDvqwb9xpsHv3m/Y58PW/79q+g6u4/MoB9S5TUk4YqiVJuTR8OMz9rUFcwQ8YyzZuvfAnvOnZ448+yoB7ftg1XHbmJ5g6tQGFSsoFQ7UkKbc+9jH4ORfxAqO4/faS1R/t7aTrr+/a719ZyG9+8vTGFCkpFwzVkqTcuvxyur54uHEjLF4M7N0L738/8dBDABzkBP6s3y38/u83rk5J2WeoliTl1qBB8Ad/UBzfcAO0P7qc9OCDXa99mU8z/+qpjBjRgAIl5YahWpKUa5/5DAwc2AHA6yvXcOCiS4nCM8m/yGf4qwGf58YbG1mhpDzo1+gCJEmqpVNPhQ9/+Bm+/e3JbGMs/+nAbQznFZYyl2Wcw2dvgDFjGl2lpKwzVEuScu/3fu9ZNm+ezOLFTdzBH3W9fsUV8PnPN7AwSbnh8g9JUu7175+4/3646SYYPRre+c7O7e9+F/r2bXR1kvLAK9WSpONCv37wuc91/khStXmlWpIkSaqQoVqSJEmqkKFakiRJqpChWpIkSapQQ0J1RNwUEdsiYnnh55JG1CFJkiRVQyPv/nFbSumWBh5fkiRJqgqXf0iSJEkVamSovjYiVkbEP0TESQ2sQ5IkSapIpJRqM3HEImBUD2/dCDwMvAQk4AvA6JTSxw8zz9XA1QAjR448+6677qpJvcqmtrY2hg4d2ugy1MvYF+rOnlB39oR6UtoX8+fP/3VKaU65n61ZqC67gIiJwL0ppTOOtm9LS0tat25dzWtSdrS2tjJv3rxGl6Fexr5Qd/aEurMn1JPSvoiIYwrVDfmiYkSMTik9VxheDqwu53NPPfVUW0SYqlVqBJ3/6iGVsi/UnT2h7uwJ9aS0LyYcywcbdfePr0TELDqXf2wGPlnm59Ydy98YlH8R8ag9oe7sC3VnT6g7e0I9qaQvGhKqU0ofacRxJUmSpFrwlnqSJElShbIWqr/Z6ALU69gT6ol9oe7sCXVnT6gnb7svGn73D0mSJCnrsnalWpIkSep1emWojoj3R8S6iNgQEZ/p4f0BEfEvhfd/VbjXtXKsjJ64ICIei4j2iLiiETWqvsroif8WEWsKT25dHBHHdGskZVMZfXFNRKyKiOURsSQiZjaiTtXP0XqiZL8rIiJFhHcEOQ6Uca74aETsKJwrlkfEHx9tzl4XqiOiL/B14GJgJvChHk56nwB2pZROA24DvlzfKlVPZfbEFuCjwPfrW50aocyeeByYk1I6C/gB8JX6Vql6K7Mvvp9SOjOlNIvOnri1zmWqjsrsCSJiGHAd8Kv6VqhGKLcvgH9JKc0q/HzraPP2ulANnANsSCltSikdBO4EPtBtnw8A3yls/wC4MCKijjWqvo7aEymlzSmllcChRhSouiunJ/4tpbS3MHwYGFfnGlV/5fTFnpLhEDqfl6D8KidTAHyBzr9k7a9ncWqYcvvimPTGUD0WeLZkvLXwWo/7pJTagd1Ac12qUyOU0xM6vhxrT3wC+GlNK1JvUFZfRMSnImIjnSHqujrVpsY4ak9ExGxgfErp3noWpoYq93fIBwtLCH8QEeOPNmlvDNU9XXHufiWhnH2UH/7/Vndl90RE/CEwB7i5phWpNyirL1JKX08pTQE+DfzPmlelRjpiT0REHzqXkd5Qt4rUG5RzrrgHmFhYQriI4gqJw+qNoXorUPq3gXHA9sPtExH9gCbg5bpUp0Yopyd0fCmrJyJiAXAjsDCldKBOtalxjvVccSfwOzWtSI12tJ4YBpwBtEbEZuBc4G6/rJh7Rz1XpJR2lvze+Hvg7KNN2htD9TJgakRMioj+wJXA3d32uRv4o8L2FcADyRtu51k5PaHjy1F7ovBPut+gM1C/2IAaVX/l9MXUkuFvA+vrWJ/q74g9kVLanVIakVKamFKaSOf3LxamlB5tTLmqk3LOFaNLhguBJ482ab+qllgFKaX2iLgW+H9AX+AfUkpPRMTngUdTSncD3wb+KSI20HmF+srGVaxaK6cnIuJdwP8FTgIui4g/Tymd3sCyVUNlniduBoYC/6fwPeYtKaWFDStaNVdmX1xb+BeM14FdFC/QKIfK7AkdZ8rsi+siYiHQTmfW/OjR5vWJipIkSVKFeuPyD0mSJClTDNWSJElShQzVkiRJUoUM1ZIkSVKFDNWSJElShQzVkiRJUoUM1ZLUi0REc0QsL/w8HxHbSsYP1eiYsyPiW0d4f2RE3F+LY0tSXvS6h79I0vEspbQTmAUQETcBbSmlW2p82M8Cf3GEmnZExHMR8Z6U0oM1rkWSMskr1ZKUERHRVvjvvIj4RUTcFRFPRcSXIuLDEfFIRKyKiCmF/UZGxA8jYlnh5z09zDkMOCultKIwfm/JlfHHC+8D/Bj4cJ3+qJKUOYZqScqm3wCuB84EPgJMSymdA3wL+M+Ffb4G3JZSehfwwcJ73c0BVpeM/xT4VEppFnA+sK/w+qOFsSSpBy7/kKRsWpZSeg4gIjYCPyu8vgqYX9heAMyMiDc+c2JEDEspvVoyz2hgR8n4QeDWiPge8KOU0tbC6y8CY6r/x5CkfDBUS1I2HSjZPlQyPkTx3N4HmJtS2sfh7QMGvjFIKX0pIu4DLgEejogFKaW1hX2ONI8kHddc/iFJ+fUz4No3BhExq4d9ngROK9lnSkppVUrpy3Qu+ZheeGsab14mIkkqYaiWpPy6DpgTESsjYg1wTfcdClehm0q+kPhfImJ1RKyg88r0Twuvzwfuq0fRkpRFkVJqdA2SpAaKiP8KvJpSOtK9qn8JfCCltKt+lUlSdnilWpL0t7x5jfabRMRI4FYDtSQdnleqJUmSpAp5pVqSJEmqkKFakiRJqpChWpIkSaqQoVqSJEmqkKFakiRJqtD/B0yeButlWIp9AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 864x360 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "CPU times: user 1.42 s, sys: 20 ms, total: 1.44 s\n",
      "Wall time: 833 ms\n"
     ]
    }
   ],
   "source": [
    "%%time\n",
    "\n",
    "dx   = 1.0   # grid point distance in x-direction (m)\n",
    "dz   = dx    # grid point distance in z-direction (m)\n",
    "\n",
    "dt = 0.001   # time step (s)\n",
    "\n",
    "mux, muz = FD_2D_SH_JIT(dt,dx,dz,f0,xsrc,zsrc)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## What we learned:\n",
    "\n",
    "* Implementation of the 2D isotropic elastic SH problem based on our 2D acoustic FD code \n",
    "* Comparison with analytical solution for homogeneous Vs model"
   ]
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.6.5"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}
